function [] = ofdm_basic_cho
% modified from "OFDM_basic.m" in 
% "MIMO-OFDM Wireless Communications with MATLAB", Yong Soo Cho
%
% ####################################################################
% following function is changed to local:
% "guard_interval.m", "remove_GI.m", "plot_ber.m", "ber_QAM", "Q.m"
% ####################################################################
%

% NgType = 1/2 for cyclic prefix/zero padding
NgType = 1;
if NgType == 1
    nt = 'CP';
elseif NgType == 2
    nt = 'ZP';
end

% Ch = 0/1 for AWGN/multipath channel
Ch = 1;
if Ch == 0
    chType = 'AWGN'; Target_neb = 100;
else
    chType = 'CH'; Target_neb = 500;
end

% figure(Ch + 1);

% Channel tap power profile 'dB'
PowerdB = [0 -8 -17 -21 -25];

% Channel delay 'sample'
Delay = [0 3 5 6 8];

% Channel tap power profile 'linear scale'
Power = 10 .^ (PowerdB / 10);

% Chanel tap number
Ntap = length(PowerdB);

% Channel length
Lch = Delay(end) + 1;

% Modulation order = 2/4/6 for QPSK/16QAM/64QAM
Nbps = 4;
M = 2 ^ Nbps;

% FFT size
Nfft = 64;

% Guard interval length (Ng = 0 for no GI)
Ng = 3;
% Ng = Nfft/4;

% Symbol duration
Nsym = Nfft + Ng;

% Nvc = 0: no virtual carrier
Nvc = Nfft / 4;
Nused = Nfft - Nvc;

EbN0 = 0 : 5 : 20; % EbN0

% Number of iterations for each EbN0
N_iter = 1e5;

% Number of symbols per frame
Nframe = 3;

% Signal power initialization
sigPow = 0;

file_name = ['OFDM_BER_' chType '_' nt '_' 'GL' num2str(Ng) '.dat'];
fid = fopen(file_name, 'w+');

norms = [1 sqrt(2) 0 sqrt(10) 0 sqrt(42)]; % BPSK 4-QAM 16-QAM

for i = 0 : length(EbN0)
    
%     randn('state', 0); rand('state', 0);
    
    % BER initialization
    Ber = 0; % Ber2 = ber();
    
    % Initialize the number of error/total bits
    Neb = 0;
    Ntb = 0;
    
    for m = 1 : N_iter
        
        % Tx
        X = randi(M - 1, 1, Nused * Nframe); % bit: integer vector
        %       X = randint(1, Nused * Nframe, M); % bit: integer vector
        
        Xmod = qammod(X, M, 'gray') / norm(Nbps);
        %       Xmod = qammod(X, M, 0, 'gray') / norms(Nbps);
        
        if NgType ~= 2
            x_GI = zeros(1, Nframe * Nsym);
        elseif NgType == 2
            x_GI = zeros(1, Nframe * Nsym + Ng);
            % Extend an OFDM symbol by Ng zeros
        end
        
        kk1 = 1 : Nused / 2; 
        kk2 = Nused / 2 + 1 : Nused; 
        kk3 = 1 : Nfft; 
        kk4 = 1 : Nsym;
        
        for k = 1 : Nframe
            if Nvc ~= 0
                X_shift= [0 Xmod(kk2) zeros(1, Nvc - 1) Xmod(kk1)];
            else
                X_shift= [Xmod(kk2) Xmod(kk1)];
            end
            x = ifft(X_shift);
            x_GI(kk4) = guard_interval(Ng, Nfft, NgType, x);
            kk1 = kk1 + Nused; kk2 = kk2 + Nused; kk3 = kk3 + Nfft; kk4 = kk4 + Nsym;
        end
        
        if Ch == 0
            % No channel
            y = x_GI;  
        else
            % Multipath fading channel
            channel = (randn(1, Ntap) + 1i * randn(1, Ntap)) .* sqrt(Power / 2);
            h = zeros(1, Lch); h(Delay + 1) = channel; % cir: channel impulse response
            y = conv(x_GI, h);
        end
        
        if i == 0 % Only to measure the signal power for adding AWGN noise
            y1 = y(1 : Nframe * Nsym); sigPow = sigPow + y1 * y1';
            continue;
        end
        
        % Add AWGN noise
        snr = EbN0(i) + 10 * log10(Nbps * (Nused / Nfft)); % SNR vs. Eb/N0 by Eq.(4.28)
        noise_mag = sqrt((10 .^ (-snr / 10)) * sigPow / 2); % N0 = Eb / SNR
        y_GI = y + noise_mag * (randn(size(y)) + 1i * randn(size(y)));
        
        % Rx
        kk1 = (NgType == 2) * Ng + (1 : Nsym); kk2 = 1 : Nfft;
        kk3 = 1 : Nused; kk4 = Nused / 2 + Nvc + 1 : Nfft; kk5 = (Nvc ~= 0) + (1 : Nused / 2);
        if Ch == 1
            H = fft([h zeros(1, Nfft - Lch)]); % Channel frequency response
            H_shift(kk3) = [H(kk4) H(kk5)];
        end
        
        for k = 1 : Nframe
            Y(kk2) = fft(remove_GI(Ng, Nsym, NgType, y_GI(kk1)));
            Y_shift = [Y(kk4) Y(kk5)];
            if Ch == 0
                Xmod_r(kk3) = Y_shift;
            else
                Xmod_r(kk3) = Y_shift ./ H_shift;  % Equalizer - channel compensation
            end
            kk1 = kk1 + Nsym; kk2 = kk2 + Nfft; kk3 = kk3 + Nused; kk4 = kk4 + Nfft; kk5 = kk5 + Nfft;
        end
        
        X_r = qamdemod(Xmod_r * norm(Nbps), M, 'gray');
%         X_r = qamdemod(Xmod_r * norms(Nbps), M, 0, 'gray');
        
        Neb = Neb + sum(sum(de2bi(X_r, Nbps) ~= de2bi(X, Nbps)));
        Ntb = Ntb + Nused * Nframe * Nbps;  % [Ber, Neb, Ntb] = ber(bit_Rx, bit, Nbps);
        if Neb > Target_neb 
            break; 
        end
    end
    
    if i == 0
        sigPow = sigPow / Nsym / Nframe / N_iter;
        fprintf('Signal power = %11.3e\n', sigPow);
        fprintf(fid, '%%Signal power = %11.3e\n%%EbN0[dB]       BER\n', sigPow);
    else
        Ber = Neb / Ntb;
        fprintf('EbN0 = %3d[dB], BER = %4d / %8d = %11.3e\n', EbN0(i), Neb, Ntb, Ber);
        fprintf(fid, '%d\t%11.3e\n', EbN0(i), Ber);
        if Ber < 1e-6
            break;  
        end
    end
end

if fid ~= 0
    fclose(fid);
end

disp('Simulation is finished');

plot_ber(file_name, Nbps);

end

%%
function y = guard_interval(Ng, Nfft, NgType, ofdmSym)

if NgType == 1
    y = [ofdmSym(Nfft - Ng + 1 : Nfft) ofdmSym(1 : Nfft)];
elseif NgType == 2
    y = [zeros(1, Ng) ofdmSym(1 : Nfft)];
end

end

%%
function y = remove_GI(Ng, Lsym, NgType, ofdmSym)

if Ng ~= 0
    if NgType == 1  % cyclic prefix
        y = ofdmSym(Ng + 1 : Lsym);
    elseif NgType == 2 % cyclic suffix
        y = ofdmSym(1 : Lsym - Ng) + [ofdmSym(Lsym - Ng + 1 : Lsym) zeros(1, Lsym - 2 * Ng)];
    end
else
    y = ofdmSym;
end

end

%%
function plot_ber(file_name, Nbps)

figure;

EbN0dB = 0 : 1 : 30;
M = 2 ^ Nbps;
ber_AWGN = ber_QAM(EbN0dB, M, 'AWGN');
ber_Rayleigh = ber_QAM(EbN0dB, M, 'Rayleigh');
semilogy(EbN0dB, ber_AWGN, 'r:');
hold on;
semilogy(EbN0dB, ber_Rayleigh, 'r-');
a = load(file_name);  
semilogy(a(:, 1), a(:, 2), 'b--s');  
grid on;
legend('AWGN analytic','Rayleigh fading analytic', 'Simulation');
xlabel('EbN0[dB]');
ylabel('BER'); 
axis([a(1, 1) a(end, 1) 1e-5 1]);

end

%%
function ber = ber_QAM(EbN0dB, M, AWGN_or_Rayleigh)
% Find ananlytical BER of Mary QAM in AWGN or Rayleigh channel
%
% EbN0dB = EbN0dB: Energy per bit-to-noise power[dB] for AWGN channel
%        = rdB: Average SNR(2 * sigma Eb/N0)[dB] for Rayleigh channel
% M = Modulation order (Constellation size)

N = length(EbN0dB);
sqM = sqrt(M);
a = 2 * (1 - power(sqM, -1)) / log2(sqM);
b = 6 * log2(sqM) / (M - 1);

if nargin < 3
    AWGN_or_Rayleigh = 'AWGN';
end

if lower(AWGN_or_Rayleigh(1)) == 'a'
    ber = a * Q(sqrt(b * 10 .^ (EbN0dB / 10)));
else
    rn = b * 10 .^ (EbN0dB / 10) / 2; 
    ber = 0.5 * a * (1 - sqrt(rn ./ (rn + 1)));
end

end

%%
function y = Q(x)

% erfc: complementary error function
y = erfc(x / sqrt(2)) / 2;

end




